Critical exponents from cluster coefficients 
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For a large class of repulsive interaction models, the Mayer cluster integrals can be transformed 
into a tridiagonal real symmetric matrix R mn , whose elements converge to two constants. This 
allows for an effective extrapolation of the equation of state for these models. Due to a nearby 
(nonphysical) singularity on the negative real z axis, standard methods (e.g. Pade approximants 
based on the cluster integrals expansion) fail to capture the behavior of these models near the 
ordering transition, and, in particular, do not detect the critical point. A recent work (Eisenberg 
and Baram, PNAS 104, 5755 (2007)) has shown that the critical exponents a and a' , characterizing 
the singularity of the density as a function of the activity, can be exactly calculated if the decay of the 
R matrix elements to their asymptotic constant follows a 1/n 2 law. Here we employ renormalization 
arguments to extend this result and analyze cases for which the asymptotic approach of the R matrix 
elements towards their limiting value is of a more general form. The relevant asymptotic correction 
terms (in RG sense) are identified and we then provide a corrected exact formula for the critical 
exponents. We identify the limits of usage of the formula, and demonstrate one physical model 
which is beyond its range of validity. The new formula is validated numerically and then applied to 
analyze a number of concrete physical models. 



I. INTRODUCTION 

It is often said that the mechanism underlying phase 
transitions is the decrease of internal energy in the or- 
dered phase. However, it has been shown long ago that 
melting is dominated by the strong short ranged repulsive 
forces, and the related solid-fluid transitions are entropy- 
driven. Accordingly, purely repulsive models have been 
often used to study the fluid equation of state towards the 
structural ordering transition. The most striking demon- 
stration of these observations is given by the family of 
hard-core models, which have long played a central role 
in this field. In these models, particles interact exclu- 
sively through an extended hard core, and there is no 
temperature scale associated with the potential (interac- 
tion energy is cither infinite inside the exclusion region 
or zero outside). Thus, temperature and energy play no 
role, and the dynamics is completely determined by en- 
tropy considerations. Yet, these models exhibit various 
types of ordering transitions. They include, for exam- 
ple, the famous isotropic-nematic transition in a three 
dimensional system of thin hard rods [H, H, as well as 
the extensively studied hard spheres models]!], 0, [H, @, 0] , 
undergoing a first order fluid-solid transition for d > 3 
and, presumably, a second order transition from a fluid 
to the hexatic phase [j| [§]. These models are purely 
entropy-driven, yet they capture the essential molecular 
mechanism that drives freezing transitions. 

A complete description of the fluid phase is provided 
by the Mayer cluster series in terms of the activity, 
z = exp(/3/i), where /i is the chemical potential. For 
purely repulsive potentials, the radius of convergence of 
the cluster series is known to be determined by a singu- 
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larity on the negative real axis, z = — zq, typically very 
close to the origin Near this point, the singular part 
of the density is characterized by the critical exponent a: 

Psing(z) ~ (z + Z ) a . 

As a result of this singularity, the radius of convergence 
of the Mayer series includes only the extremely low den- 
sity regime, and the fluid-solid transition is way beyond 
it. It is therefore desirable to find a way to extend the 
information contained in the cluster integrals series to 
provide information about the behavior of the system 
close to the ordering transition region. In particular, one 
is interested in the critical exponent a' characterizing the 
density near the physical termination point of the fluid 

Psing(z) ~(z- Zt)" 7 . 

It has been shown that this goal may be achieved by 
transforming the cluster integral series into a tridiago- 
nal symmetric matrix form The matrix elements 
Rnm adopt a clear asymptotic form, and converge ex- 
tremely fast to two different constants: A (off-diagonal) 
and B (diagonal). This fact can then be utilized to obtain 
good approximants for the fluid density far outside of the 
convergence circle of the power series [S mEil. Like 
Pade methods, these approximants are consistent with 
the known elements to all available orders. However, the 
R matrix scheme seems to fit much better purely repul- 
sive systems, as it incorporates the existence of two sin- 
gular points on the real axis [HI, [l7|. Yet, a major 
shortcoming of this approach was its failure at the critical 
regime. It is easy to prove (see below) that tridiagonal 
R matrices described at the asymptote by two constant 
values lead to universal critical exponents a = a' = 1/2 
at both singularities, which are obviously wrong. Thus, 
the above approach fails when one is in close vicinity to 
the transition region. 
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TABLE I: Mayer cluster coefficients nb n for various models 



A partial solution for this problem was recently found, 
noticing that for many of the studied models not only the 
matrix element approach a constant but also the asymp- 
totic correction to the constant takes a universal form, 
following a 1 /n 2 decay of the elements to their constant 
asymptotic value [l8[: 

B„ = Ra,„ = B + b/n 2 
A n+1/2 = Rn, n+ i = A + a/(n+ 1/2) 2 . (1) 

Under these circumstances, one is able to analytically 
calculate the critical exponents at both fluid termina- 
tion point (the physical one, at the ordering transition 
or at the termination of the super-cooled fluid, and the 
nonphysical one on the negative real z-axis). These expo- 
nents depend on the amplitudes of the 1/n 2 corrections, 
and generally deviate from 1/2. This approach works 
satisfactorily for many models and tests well against the 
known result for the nonphysical singularity that predicts 
universal critical exponents depending on dimensionality 
alone. Yet, while many models indeed show this simple 
1/n 2 decay, we have found out that some other models 
exhibit different asymptotic behavior. For exam ple, the 
R matrix elements of the hard hexagons model are 
presented in fig [TJ As this is an exactly solvable model, 
one is able to produce a large number of cluster inte- 
grals. Doing so, we note that while the first few elements 
seem to follow the 1/n 2 rule, the asymptotic behavior is 
quite different. The matrix elements do converge to two 
constants as expected, but their leading asymptotic be- 
havior follows an oscillatory 1/n decay rather than the 
above mentioned 1/n 2 . This finding raises the question 
of how to deal with R matrices whose correction deviates 
from the |T]) form. Moreover, it sheds doubt on the ap- 
plicability of former results to other models where only 
a few cluster integrals are known: one may argue that 
the hard hexagons example shows that the 1/n 2 behav- 



ior is only a transient one, and the true asymptotics of all 
these models is different. Indeed, extension of the avail- 
able series to higher coefficients of the Mayer expansion 
allowed us to see in a number of additional models that 
the seeming 1 /n 2 behavior is accompanied by additional 
corrections, including an oscillatory cos(qn)/n term that 
becomes dominant in the asymptote. We observed such 
oscillations, for example, for hard-core two-dimensional 
square lattice gas with exclusion shell up to second (N2 
model) , third (N3 model) and fourth (N4 model) nearest 
neighbors. 

As this oscillatory term dominates for large n, the va- 
lidity of the results of [18| is put in question. Therefore, 
we set out to study the effect of this additional correc- 
tion term on the critical behavior of the equation of state. 
Here we extend the previous result and explore the case 
of matrix elements taking the asymptotic form 



Rn,n 
Rn^n+l 



B+ b +b ,co S ( q n) 
n n 



A 



a 



(n + 1/2) 2 



, cos[ g (n+l/2)] 
a n + 1/2 (2) 



Using an analytical RG-like decimation scheme, we show 
that in this case the critical exponents are given by 



1 / A(2a + b) 
a= 2\ l — 



[2a' cos(<j/2) + b 1 ] 2 
[1 - cos( g )]A 2 



, 1 4(2a - b) 

° =2\ 1 -^— 



[2a 1 cos(g/2) - b'} 2 
[1 - cos(q)]A 2 



(3) 



thus generalizing the results of [18|. We also discuss the 
possible effect of other kinds of corrections, and conclude 
that they do not affect the critical exponent as long as 
the spectrum of the matrix remains intact. We verify 
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(4) 



where b n is the n th Mayer cluster integral. It is always 
possible (see appendix A for an explicit construction) to 
define a tridiagonal symmetric R matrix which satisfies 
the condition (n > 1) 



(R n )n = (-l) n (n + l)b n+1 . 
The density may then be expressed in terms of R: 
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FIG. 1: Hard hexagons matrix elements, fitted to the form 
©, with A = %/l25/4 a = 0.0027 a' = -0.063 B = 5.5 
b = 0.627 b' = 0.129 q = 0.36 



the result by extensive numerical study of artificial mod- 
els and by analysis of the exactly solvable hard-hexagons 
model. The next-nearest neighbor exclusion model on a 
triangular lattice is discussed as an example in which the 
spectrum does not remain intact and our approach breaks 
down. Finally, we apply our formula to the two models 
that have been recently studied by means of Monte-Carlo 
(MC) simulations [20(: the hard-core two-dimensional 
square lattice gas with exclusion shell up to fourth (N4 
model) and fifth (N5 model) nearest neighbors. 



II. ANALYSIS 

For the sake of completeness, we start with a brief 
review of the approach presented in [l8l |. The Mayer 
cluster integrals provide a low-z expansion for the density 



p(z) = nb n z n = J2 {-l) n z n+1 {R n )n = z{I + zR)^. 



n=l 



n=0 



(6) 

Alternatively, the matrix inversion in the previous equa- 
tion may be expressed in terms of the spectrum A of the 
R matrix, and the corresponding eigenvectors tp(X): 



^i(A) 2 



z- 1 +A' 



(7) 



where ipi(X) is the first component of the ip(X) vector. 

The reciprocals of the eigenvalues of this matrix are the 
Yang-Lee zeroes of the grand-canonical partition func- 
tion. For all purely repulsive models studied to date, the 
R matrices are real- valued, and thus their eigenvalues are 
also real (R is symmetric by construction). There is yet 
no proof that this is indeed the case for all such mod- 
els, but construction of R matrices for dozens of different 
lattice and continuum purely repulsive models (see, e.g., 
0) El) Ell and this work) provides strong evidence for 
it: in all cases studied the matrix elements were real to 
all orders calculated. Furthermore, as mentioned above, 
the matrix elements in all models studied adopt a clear 
asymptotic pattern, converging quickly to a (real) con- 
stant. Therefore the possibility that some higher order 
element may become complex seems improbable. 

For these real R matrices the spectrum of the matrix 
lies on the real axis in an interval (— z t _1 , Zq 1 ) (and the 
Yang-Lee zeroes lie on two intervals along the real activ- 
ity axis: z < — zq and z > z t ). It follows from ([7]) that 
the density p(z) has two singular points at z values for 
which — z^ 1 coincides with the spectrum edges of the R 
matrix, leading to vanishing of the denominator on the 
right-hand side. The critical behavior of the density p(z) 
near the physical and non-physical singularities is there- 
fore determined by the structure of the residue ^i(A) at 
the spectrum edges. 

For example, we look at a matrix with two con- 
stants along the three main diagonals, B (diagonal) 
and A (off-diagonal). The eigenvalues are A(fc) = 
B + 2Acos(k) (0 < fc < 7r) and the eigenvectors are 
ijj n [\(k)] = sin(nfc). The critical points are then 

-z^ 1 = -{B + 2A) 



4 



(corresponding to k = 0), and 



= 2A-B 



where cr(cr') is the exponent of the non-physical (physical) 
branch point. 



[k = 7r), where ipi(k) = ipi[X(k)] ~ k and ipi(k) = 
ipi[X(k)] ~ (fc — 7r) respectively. Expanding the integral 
in ([7]) for z ~ — zq and z ~ zt one finds that the density 
terminates at both ends with a square-root singularity. 
We now consider a general R matrix taking the form 



B n = R n ,n — B + 8B n 
A n +l/2 = Rn,n+l = A + 5A n+ i/ 2 - 



(8) 



The critical behavior is determined by the long- 
wavelength, slowly-varying, eigenvectors and therefore 
the eigenvalue equation (we treat the nonphysical crit- 
ical point only, analysis of physical point is essentially 
identical) 

A n -l/2lpn-l + B n 1p n + A n+1 / 2 1p n +l = Xlp n (9) 

may be studied in the continuum limit, taking the form 
of a differential equation in the variable x = kn. For the 
general case , the discrete equations © transform into 



A*) + /(*) + [ " n -^\ " +1/2J m = o. 

(10) 

As long as the corrections SB and SA are small enough 
(see below) the spectrum does not change. The eigenvec- 
tors, nevertheless, are modified. In [l8[ the R matrix was 
assumed to take the form (fTJ, and then the differential 
equation (|10|l is reduced into a Bessel equation. A closed 
form for the eigenvectors is available, and one obtains 
the critical behavior of the density near the two branch 
points p(z c ) — p{z) = (z c — z) a (or p(z) = [z c — z)~ a if the 
density diverges at criticality, such as the case of the non- 
physical singularity in d < 2). The critical exponents are 
given by 



2a - b 
A ' 



^2a + b 
A 



This approach, however, cannot be extended straight- 
forwardly to study a general correction to the matrix 
elements: while for 1/n 2 corrections (flQ|) can be written 
in terms of x = kn alone, independently of k, a general 
correction term results in a fc-dependent differential equa- 
tion. More importantly, considering terms 0(l/n 3 ) in the 
differential equation approach leads to an essential singu- 
larity at the origin, resulting in transition layer solutions 
and complicated behavior at the origin. These terms in- 
deed show up when one analyzes real R matrices (see 
below for the N4 and N5 models). Third, the mapping 
to a differential equation relies on the slow variation of 
the eigenvectors and is bound to fail for correction terms 
of the form ([2]) that induce an intrinsic "length" -scale (on 
the n axis) into the problem. 

We thus present here a complementary approach to 
study the general correction term, which is based on the 
idea of renormalization. In their discrete form, the eigen- 
value equations ([9]) form an infinite linear system of equa- 
tions. Since the system is tridiagonal, it is quite easy to 
eliminate half of the variables, e.g. all variables ip n for n 
even. This effectively removes half of the rows and half 
of the columns in the matrix, "tracing out" half of the 
degrees of freedom in the problem. One obtains a new 
tridiagonal system of equations, or a renormalized R ma- 
trix, with the same eigenvalues and new vectors ip(k) that 
are simply related to the former ones ip n (k) — ip2n-i(k). 
In particular, ipi(k) = ip\{k). The density as a func- 
tion of z is fully determined by the spectrum and ipi(X) 
through ([7])- Thus, the renormalized R matrix may be 
utilized to generate the same equation of state and the 
same critical behavior as the original one. Explicitly, the 
reduced eigenvalue equation after one such decimation 



(11) 



process takes the form (n odd; A Ti 
n<0) 



-1/2 



= B n = for 



X — B n _i 



1pn-2 + ( 



A 2 

X — B n _i 



A 2 

X — B n+ \ 

r 



B n )tpr, 



^n+1/2 Ai+3/2 

X — B n+ \ 



tpn+2 = Xljjj, 



(12) 



Accordingly, the R matrix elements transform, under 
such decimation, according to 



BL 



A 2 

^2ra-3/2 



.4 



n+1/2 



A 2 

^2n-l/2 

X — B<z n -2 X — B-2n 
^-2ri-l/2^-2n+l/2 



B 2n -i (13) 
(14) 



X — B 2n 

In the transformed linear system %j> n is in fact V^n-i, so 



for a given functional form for A n and B n one should 
change variables n' — > 2n — 1. Note that the renormal- 
ization transformation is A-dependent. Since the density 
in the vicinity of the critical points is determined by the 
spectrum edges only, this poses no difficulty. 

As a first demonstration of this RG scheme, one may 
look at the solvable case of 1 Jn 2 correction. Substituting 
A n = A + a/n 2 , B n = B + b/n 2 and X = 2A + B into 



5 



([12]) . one obtains 

A n+l/2 - A/2+-L(4a + 6) + 0(i) ( 15 ) 

' Sn z n J 

B„ -► (A + S) + ^-(a/2 + 36/8) + 0(^r). 

n z vr 

Clearly, the spectrum edge, defined by the asymptotic 
value of A n _x/ 2 + A n+ x/ 2 + B n to be — zq = — (2A + B)^ 1 
is conserved under the decimation. Moreover, the correc- 
tion term (5A n _i/2 + 5A n+ i/2+SB n )/A which appears in 
the differential equation and determines the critical ex- 
ponent by pip , is also stable under the transformation 
and remains equal to (2a + 6) /An 2 , as expected. 

Applying the same transformation for corrections of 
the form l/n Q i.e. A n = A + a"/n a ,B n = B + b"/n a , 
results in 

1 , r „ S A ^ 1 2a" + 6" 1 
■j{SA n _ 1/2 + SA n+1/2 + SB n ) -» — — . 

(16) 



Therefore, one may conclude that for a > 2 the correc- 
tion term in the differential equation IjlUp is suppressed 
by successive applications of the RG decimation transfor- 
mations. Therefore, these correction terms are irrelevant 
in determining the critical exponents. 

We now employ the RG scheme to study the case of 
main interest: 1/n- modulated oscillations, as observed 
for the hard hexagons model 



A n+1/2 = A + a'cos[g(n + l/2)]/(n+l/2) 

B n = B + b' cos(gn)/n. (17) 



The transformation of the differential equation correction 
term (SA n _ 1 / 2 + dA n+1 / 2 + SB n )/A upon one decimation 
step is given by 



-j(SA n _y 2 + 8A n+1/2 + 5B n ) 



2a! cos(q/2) + b' 
A 



[2a'cos(g/2) + b'][l + cos(g)] cos(2gn) [2a' cos(g/2) + b'] 2 1 



An 



2A 2 



+ 0(cos(2g)/n 2 ,l/n 3 



(18) 



Obviously, the real-space renormalization process in- 
duces a change in the frequencies q — > 2g. In addi- 
tion, (i) the cos(gn)/n term is multiplied by a factor of 
[1 + cos(g)], and three more terms emerge: (ii) a new 
1/n 2 term, (iii) terms 0(cos(2gn)/n 2 ) and (iv) terms 
0(l/n 3 ). Iterating this procedure N times, one obtains 

N 

from (i) (2a + 6) -> (2a' + 6') Y\[l + cos(2 n - 1 q)}. The 

n=l 

newly emerging 1/n 2 terms (ii) combine to take the form 
(2 °+j } J2 II [l+cos(2 m " 1 g)] 2 . The first term gets ex- 

n— m— 1 

N 

ponentially small for large N: Y[ [1 + cos(2" _1 g)] ~ 4~ N 

(see Appendix B), and thus could be neglected. The 
sum over the products in the second term converges to 
[1 — cos(g)] -1 (see Appendix B). This second term does 
affects the critical behavior as it adds up to the l/?i 2 
terms in the R matrix. The 1/n 3 terms (iii) may be ne- 
glected as their amplitude decreases: each existing 1/n 3 
term decreases by factor 2 upon an RG step, accord- 
ing to (|16p. While a new term is being added from the 
transformation (fT2|) , the sum of all contributions still de- 
creases exponentially with the number TVof RG steps. 
The cos(gn)/n 2 terms (iv) transform under decimation 
in an analogous way to the original cos(qn)/n term: they 
get multiplied by a factor 1 + cos(g) resulting in an expo- 



nential decay, and give rise to new 0(1/ n ) terms (analo- 
gous to the 0(l/n 2 ) terms generated from decimation of 
the cos(qn)/n), as well as faster decreasing terms. Again, 
the 1/n 4 exponentially decrease through decimation by 
(fTBf and are therefore neglected. In summary, the net 
effect of the cos(qn)/n term after a large number of RG 
steps is the creation of a new 1/n 2 term. These terms, 
emerging from the decimation process, can then be ana- 
lyzed using the mapping to the Bessel differential equa- 
tion as described in |l8|. 

Up to this point we treated the pure cos(gn)/n case. 
Similar analysis may be done for the mixed case, where 
both cos(qn)/n and 1/n 2 terms are present (as happens 
for the physical models to be discussed). It turns out that 
the transformation equation (|12p does mix the correction 
terms, as the numerator of [A n A m )/(\ — B{) is quadratic 
in the off-diagonal elements. However, the mixed terms 
will follow the form cos(gn)/n 3 which can be ignored 
based on arguments similar to those presented above for 
the cos(gn)/n 2 terms. Multiplicative cross-terms could 
be relevant (in RG sense) only if they decay 0(l/n 2 ) 
or slower. Thus, one may simply add the original 1/n 2 
terms to those emerging from RG. Collecting the 0(l/n 2 ) 
terms originating from both the functional form of the 
matrix elements and the decimation process for oscilla- 
tory terms, one obtains the closed form §3§ for the criti- 
cal exponents in the general case (having both cos(gn)/n 
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FIG. 2: Critical exponents as measured from equation of 
state (|19|) . compared with the exact prediction ([3j> for var- 
ious choices of cos(qn)/n corrections. 



Finally, we note that if the amplitude of the correc- 
tion to a constant matrix is strong enough, one obtains 
from ([3]) an imaginary value for a. When this happens, 
both solutions of the differential equation PH|) diverge at 
the origin [21] . Consequently, there are no solutions to 
the eigenvalue problem Q for k ~ 0, or A ~ 2A + B. In 
other words, a perturbation of the constant matrix which 
is strong enough to make a imaginary modifies the spec- 
trum of the matrix, such that the spectrum edge shifts 
from 2A + B. In these cases the critical point is not 
given by —zq = — (2 A + £?) _1 . Similarly, whenever a' 
becomes imaginary, the physical singularity shifts from 
z t = (2 A — B)^ 1 . In both cases, the corresponding crit- 
ical exponents are not given by (|3|). This scenario is re- 
alized for the next-nearest neighbor exclusion model on 
the triangular lattice (see below). 



III. NUMERICAL STUDY 



terms and 1/n 2 terms). 

We note that the critical exponents in Q do not de- 
pend on B. This can be readily understood looking at 
equation ([7]). A change in B results in a constant addi- 
tion to the whole spectrum of R, without modifying its 
eigenvectors. Looking at the expression ((TJ for the den- 
sity, it becomes clear that the effect of a constant added 
to the eigenvalues of R on the density is equivalent to a 
constant shift in z^ 1 . That is, if the density of the B = 
matrix is po(z), the density for finite B, p{z\ B), is simply 

p(z;B)=p { z _ 1 1 +B ). 

Therefore, the value of B affects the location of the crit- 
ical points, but do not change the critical exponents. 



In order to test our results, we have constructed var- 
ious tridiagonal symmetric R matrices with prescribed 
matrix-elements asymptotic form, and compared the pre- 
diction ([3]) with the critical behavior as measured from 
the equations of state calculated by ([6]) for these mod- 
els. First, we looked at matrices obeying (fT7|) . with the 
parameters A — 16 and B = 31 (z t — 1), and various 
combinations of a', b' and q. We have also modified A 
while keeping the other parameters fixed to check the A 
dependence. For these R matrices, one is able to con- 
sider as many coefficients as desired. Thus, the size of 
the sub-matrices studied is considerable, and the matrix 
inversion of ([6]) is costly. Instead, the density p can be 
equivalently calculated using the continued fraction rep- 
resentation 



Rl,l + 1/2- 



#1,2" 



R 2 ,2 + 1/z - 



■£13.3 



1/z- 



(19) 



which typically converges rather quickly (except for the 
immediate vicinity of the transition point). Figure [5] 
compares the critical exponent a 1 obtained by fitting the 
density as given by (TT9"]) for z close to the termination 
point Zt with the theoretical prediction of The re- 
sults are in excellent agreement, except for a few points 
where the numerical calculation of the density was dif- 
ficult due to slow convergence of the continued fraction 
in the immediate vicinity of the critical point. We also 
calculated the density for R matrices with both 1/n 2 and 



cos(n)/n corrections, i.e. following ([2]). The agreement 
between the theoretical prediction of ^ and the mea- 
sured critical exponent was again excellent. Another spe- 
cial case we checked was that of a 1/n 3 correction. This 
is the most dominant correction for which we predict no 
change to the critical exponent of the const matrix. Us- 
ing the same constants A and B, we looked at correction 
amplitudes up to a" — 30, and verified that the critical 
exponent indeed does not change: a' — 0.5 as expected. 
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N4 


N5 


Triangular N2 


n 


B n 




B n 


A n 


B n 




1 


21 


9.3808315196 


25 


11.489125293 


13 


6 


2 


17.045454545 


8.7248457800 


20.454545455 


10.570731201 


10.555555556 


5.5674871873 


3 


17.024781724 


8.6098449927 


20.518434015 


10.405877593 


10.550189740 


5.4798922624 


4 


17.018848337 


8.5688057487 


20.535452306 


10.348827163 


10.576974307 


5.4386602176 


5 


17.017061106 


8.5493273256 


20.540912637 


10.322722326 


10.600981921 


5.4160504788 


6 


17.016534640 


8.5385152615 


20.543142266 


10.308577639 


10.615727211 


5.4037594351 


7 


17.016426427 


8.5318795598 


20.544349327 


10.299999108 


10.622716309 


5.3973336220 


8 


17.016464632 


8.5275092336 


20.545165481 


10.294371298 


10.625196191 


5.3938592930 


9 


17.016552228 




20.545790085 









TABLE II: R matrix elements for various models 



IV. APPLICATIONS TO PHYSICAL MODELS 

Analysis of the R matrix as detailed above may be 
used to predict the critical behavior of all models with 
purely repulsive interactions. Our results apply equally 
to continuum and lattice models in all dimensions. Here 
we demonstrate applications to a number of 2D hard-core 
lattice gas models. 

For all the models to follow, we have calculated the 
cluster integrals to a high order (in order to calculate the 
R matrix). It is natural to compare standard series anal- 
ysis methods [13] to the results to be obtained from the 
i?-matrix. We have applied the ratio method, Dlog Pade 
and differential approximants to the models to follow. In 
general, ratio analysis of the series provide a rather exact 
estimate of the non-physical singularity location zq and 
the related a = 1/6, but says nothing about the physi- 
cally relevant zt and a' . Dlog Pade approximants again 
converge nicely to predict a singularity at — zq but show 
no consistent pole anywhere on the positive real z-axis. 
Similar results were obtained using the differential ap- 
proximants. Overall, these methods do better then the 
R matrix for the nonphysical singularity. The reason for 
these failures is the existence of a branch-cut singularity 
located so close to the origin, which makes the physical 



singularity, typically much further away, undetectable by 
these methods. The R matrix, which incorporates the 
branch-cut naturally, is more successful. 

Even though standard series analysis methods are of- 
ten superior to the R matrix as a means to analyze the 
non-physical singularity, we still include in the following 
the 7?-matrix results for both singularities. The reason is 
that unlike standard methods, R matrix is expected to 
work equally well for both termination points. The accu- 
racy of both exponents a and a' depends roughly equally 
on the quality of the fitting parameters describing the 
asymptotic behavior of the matrix elements. Thus, our 
R matrix results for a should not be taken as the yard- 
stick for measuring R matrix vs. Dlog Pade, but rather 
as a measure of the accuracy of the R matrix itself, as one 
expects the same degree of accuracy for both exponents 
calculated. 

A. Hard hexagons model 



The hard hexagons model (lattice gas on on a triangu- 
lar lattice with nearest-neighbors exclusion) was solved 
exactly by Baxter [l9[ . This allows us to calculate many 
cluster coefficients and matrix elements. The density in 
this model is given exactly by the relation [23| 



J 



p n (p - l)z 4 - p 5 {22p 7 - 77p 6 + 165p 5 - 220> 4 + 165p 3 - 66p 2 + 13 p - l)z 3 
p 2 {p - l) 2 (119p 8 - 476p 7 + 689p 6 - 401p 5 - 6p 4 + 125p 3 - 63p 2 + 13p - l)z 2 
(p - l) 5 (22p 7 - 77p 6 + 165p 5 - 220/ 



165p J - 66p 2 
I 



13/o- l)z + p(p- l) 11 = 0. 



(20) 



Using this relation, one is able to expand the density 
in power series of the activity z and extract the cluster 
integrals nb n . Employing infinite-precision integer com- 
putation we extended the 24 elements calculated in [23| 
to 1100 elements, enabling the construction of the first 
550 diagonal and off-diagonal elements of the R matrix 
[2~il ]. These allowed unambiguous determination of the 
asymptotic form of these elements. One can observe in 
figure [T] clear oscillations of the matrix elements. There- 



fore application of the formula presented in [18| , which is 
based on a 0(n~ 2 ) correction term, was doubtful. Based 
on the analysis above and the extended formula J3J , one 
may calculate the critical exponent from fitting the ma- 
trix elements of the hard hexagons model. This results 
in a' = 0.6662 where the exact result is a' — 2/3. Note 
that the early version of Q as presented in [l8| gives 
in this case a' = 0.6902. The result for the nonphysical 
critical exponent calculated based on our R matrix anal- 
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ysis and (J3j) is cr = 0.1655, which compares well to the 
exact universal result a = 1/6 [H, EE E3 • 

B. Triangular lattice N2 model 

Next, we study the triangular lattice N2 model (exclu- 
sion up to the next-nearest-neighbor). This model was 
long ago investigated, and early studies suggested that 
the phase transition is first order (2f| HE[l3|- However, 
later transfer matrix analysis [28| . and recent exhaustive 
MC results [29[ concluded that the model undergoes a 
second order phase transition at p c = 1.75682(2) and 
critical density p c = 0.180(4), and is believed to be part 
of the q = 4 Potts universality class, with a' = 1/3. 

We used the transfer matrix method to obtain an ex- 
act expansion of the partition function in powers of the 
activity. We have constructed transfer matrices for strips 
with width up to M — 26 (number of symmetry re- 
duced states in the M = 26 matrix is 730100). We 
then constructed the exact low-z power series expan- 
sion for the density p(z), the first 17 coefficients of which 
are identical with their bulk values (the cluster integrals 
for the models considered henceforth and the resulting 
R matrices are given in tables Q] [TT]) . The difference 
Ai-i/2 + Ai+1/2 ~ B n should converge to 2A — B — z t _1 . 
In the absence of oscillatory terms, the slope of this dif- 
ference against 1 /n 2 determines the critical exponent by 
([3]) . As seen in figure [3] the matrix elements are well fit- 
ted, with 2A - B = 0.107(1) and 2a - b = 3.61(1), and 
extrapolation of A n alone gives A = 5.382. Therefore, 
in this case analysis of the R matrix shows clearly that 
4(2a — b)/A~ 2.7 > 1 which means that (JTTJ> will lead to 
an imaginary a'. As discussed above, in such cases the 
above analysis breaks down as the spectrum edge shifts 
from 2A ± B. Indeed, for this model the critical activity 
as determined by MC studies, z t = 5.794 [2{|, deviates 
significantly from (2A—B)~ 1 = 9.35, clearly demonstrat- 
ing the spectrum edge shift. 



C. Square lattice N4 model 

Having tested the limits of the method, we move on 
to apply it and examine models in which the critical be- 
havior is not known. The N4 model on a square lat- 
tice (hard-core exclusion of all neighbors up to the 4" 1 
order) was first studied using transfer matrix methods 
[27L |3Q ] . Recently, it was revisited, employing MC sim- 
ulations [2(j. It is believed to undergo a second order 
fluid-solid transition of the Ising universality class. The 
critical chemical potential was found to be \i c = 4.705 
with a critical density p c = 0.110 [20| . where the closest 
packing density is p cp = 0.125. 

Here too, we used the transfer matrix method to ob- 
tain an exact expansion of the partition function and 
expand the density in powers of the activity. We have 
constructed transfer matrices for strips with width up to 



M = 37. Employing translational and inversion sym- 
metries, the number of symmetry reduced states in the 
M = 37 matrix is 4137859. Using this matrix, we ob- 
tained the first 18 coefficients that are identical with the 
bulk values. The diagonal matrix elements take the form 
B + b/n 2 + b' cos(qn + <f>)/n, while the off-diagonal ones 
exhibit no visible oscillations, and are well fitted by the 
cubic form A + a/n 2 + a"/n 3 (see figure @|. Based on 
the fit parameters, one is able to predict the non-physical 
singularity location — z$ = —0.0294, which compares well 
with the value we obtained from direct ratio analysis of 
the series — zq = —0.029374(1). The critical exponent at 
this singularity is calculated by[3]to be a — 0.1891, close 
to the exact universal value a = 1/6. 

Looking at the physical singularity, one observes 2 A — 
B = 0.015(5), i.e., p c = 4.2(4), barely consistent with 
the result of [2(| • While the accuracy in determining the 
critical activity is low, the critical density can be deter- 
mined to much better accuracy p c — 0.112(1), in good 
agreement with the MC results. It is remarkable that we 
are able to determine to such accuracy the critical den- 
sity at the fluid-solid transition based on the low-density 
behavior of the fluid alone. The critical exponent may 
be found by ([3]) to be a' — 0.28(6). Thus, based on our 
analysis of the cluster integrals we can quite safely ex- 
clude the possibility of the Ising universality class, where 
a' = 1 . The latter result contradicts the numerical obser- 
vations of [20j. Detailed numerical studies of this model 
aimed at an accurate calculation of the critical exponents 
are required to settle this discrepancy. 



D. Square lattice N5 model 

Finally, we look at the N5 model on a square lattice 
(hard-core exclusion of all neighbors up to the 5 th order) . 
This model was also recently studied using MC simula- 
tions [20I ] and found to undergo a weak first order transi- 
tion at fi c — 5.554. Again, we calculated 18 cluster coeffi- 
cients using the transfer matrix method up to M = 37. In 
this case, one observes no oscillations, but the R matrix 
elements exhibit a strong third-order correction term: 
A n = A + a/n 2 + a"/n 3 and B n = B + b/n 2 + b"/n 3 (see 
figure [5]). While the third order term is stronger than 
the second-order one in the regime studied, our RG anal- 
ysis allows us to conclude that the 1/n 3 correction does 
not change the critical exponents and we can use (fTT|l . 
The non-physical exponent a calculated from the above 
parameters, a — 0.1718 is in reasonable agreement with 
the exact universal result 1/6. Similar calculation for the 
physical singularity yields a — 0.1621. The accuracy of 
the latter result might suffer from the lack of insufficient 
cluster integrals. However, one can safely say that the di- 
agonal 1/n 2 amplitude b is small, and thus the physical 
exponent a' would not deviate much from cr, and should 
satisfy a' ~ 1/6. The critical activity z t = (2A — B)^ 1 
is estimated to be Zt = 166, but is highly sensitive to 
small errors in A and B and might be very well equal or 
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FIG. 3: Triangular N2 Matrix elements. The difference 
An+i/2 — B n + A n+1 / 2 extrapolates to 0.107(1), much lower 
than 1/zt = 0.1726. The slope with respect to n~ 2 is 3.61(1), 
much larger than A/4. These two observations are consistent 
with a spectrum edge shift. 



FIG. 4: N4 Matrix elements. The diagonal terms are fitted 
to the functional form ©: B = 17.0121, b = 0.19, b' = 0.029, 
q — 0.295. The off diagonal term fit well Q with an added 
cubic correction a"/n 3 : 2 A = 17.0316, 2a = 1.634, 2a = 0, 
2a" = 0.09" 



higher than the one reported in [20( (z c = 258). If z t > z c 
then the critical point we found corresponds to the ter- 
mination of the super-cooled fluid phase. This scenario 
is discussed in 13111 and was suggested to be related to a 
glass transition [l8L l3ll ] . 




O R Matrix B n 
□ R Matrix A n 
— Fitted Asymptote 



CONCLUSION 



4 6 
n 



The R matrix representation of the Mayer cluster in- 
tegrals converges very quickly to its asymptotic form. It 
therefore provides a powerful tool for extrapolating the 
low-z expansion of the fluid equation of state to cover the 
full fluid regime. In this work we analyze the analytic 
properties of this equation of state in the vicinity of the 
critical points. It is shown that not only the location of 
the critical points, but also the critical exponents can be 
determined if one identifies correctly the asymptotic be- 
havior of the R matrix elements. A number of correction 
forms are analyzed, most of which are shown by RG ar- 
guments to be irrelevant for the critical behavior. Thus, 
we provide an exact formula for the critical exponents, 
depending on a relatively few parameters characterizing 
the functional dependence of the matrix elements. Ap- 
plication of this method to a number of lattice-gas mod- 
els results in partial agreement with recent MC studies. 
Analysis of the discrepancies through an extensive MC 
study is left for future work. 



FIG. 5: N5 Matrix elements. Diagonal and off-diagonal terms 
fit well {J} with an added cubic correction b"/n 3 (a"/n 3 ) : 
B = 20.547, b = -0.0166, b" = -0.706. 2A = 20.553, 2a = 
2.28, 2a" = 0.143. 
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VI. APPENDIX A: CONSTRUCTION OF THE 
i?-MATRIX 

Here we give an explicit recursive construction of a 
tridiagonal symmetric R matrix that satisfies ([5]). First, 
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assign 



Assuming all elements Rij are known for 1 < i, j < m 
(and ([5]) is satisfied for n < 2m), we construct R m ^ m+ i 
and R m +i tTn +i as follows: 

Define P to be the m x m leading submatrix of R, i.e., 
the first to rows and first to columns of R. The next 
off-diagonal element is given by 



R 2 = 7? 2 



(2TO+l)6 2m+1 -(P 2m ) 11 



(P" 



I 2 

/ lm 



Now define Q to be the (m+1) x (to+1) leading submatrix 
of R, with zero as its m + 1, m + 1 element. The next 
diagonal element is then given by 

R -(2TO + 2)6 2m+2 -(Q 2m+1 ) 11 

W )l,m+l 

It is easy to see by explicit multiplication that the sub- 
matrix up to row and column to satisfy ([5]) up to n = 2m. 
Further matrix elements do not affect (B?)ii for j < n. 
Therefore, each additional cluster integral allows for one 
additional R matrix elements. It should be pointed out 
that the above process is exponentially sensitive to er- 
rors. This means that if one is interested in matrices 
with to > 5 or so, the cluster integrals used should be 
exact or at least known to high accuracy. In addition, 
the actual construction of R matrices should generally 
be done using high-accuracy arithmetics to avoid build- 
up of round-off errors. 



VII. APPENDIX B 



irrational, the sequence 2 J q( mod 2ir) is uniformly dense 
in (0, 2tt). Thus, in the limit N — > oo the sum may be 
replaced by an integral 



n y 

2^1n[l+cos(2 J '- 1 g)] = 27V / ln[l+cos(x)]da; = -2JVln(2). 
i=i n 



Exponentiating the result, one reveals ([21 
Secondly, we show that 



/(?) 



1 



i=i j=i 



-cos(2 j - 1 q)]' 2 



(22) 



1 — cos(q) 

It follows from the definition that f{q) satisfies /(g) — 1 = 
[1 + cos(g)] 2 /(2g). This recursion rule is indeed satisfied 
by /(g) = 2/[l — cos(g)]. All left to be shown is that there 
is no other (continuous) solution. Assume there exist 
two different solutions fi(q) and /2(g)- Their difference 
tif(<l) = fi(l) ~ h(q) then satisfies 



5f(q) = [l+cos(q)} 2 6f(2 q ) = 5f(2 n q) J] [l+cos(2^ 1 g)] 2 

3=1 

(23) 

Let g/(27r) be irrational. 5f is continuous, thus for each 
e there exists S such that \q\— q\ < 5 ^ \f(qi) — f(q)\ < e - 
Again we use the fact that the sequence 2 J q( mod 2n) is 
uniformly dense in (0, 2ir) to deduce that there exists also 
N such that \2 N q( mod 2n)-q\ < 5 and thus \6f(2 N q)- 
Sf(q)\ < e. In fact there are infinitely many such N'a, 
so one may find N as large as required to satisfy the 
latter inequality, while at the same time satisfying (|2ip . 
Employing (|23[) one finds 



We first show that 



iV 



JTjl + cos^'-^)] 2 ~4 
3=1 



— N 



(N^oo). (21) 



Taking the logarithm of the product, one obtains 

N 

2 l n [l + cos(2 J ' — 1 g)]. It is easy to see that for q/(2ir) 
3=1 



Sf(q) = 6f(2 N q)Y[[l + cos(2^ 1 9 )] 2 ~ 4T N 8f{2 N q). 
3=1 

(24) 

That is, \f{2 N q) - f(q)\ ~ (4 W - l)\f(q)\ » e in con- 
tradiction to the abode, unless Sf(q) = 0. Since this is 
true for all irrational q/(2n), the function must vanish 
identically if continuous. Q.E.D. 
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